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Abstract 

Evidence of discrete scale invariance (DSI) in daytime healthy heart 
rate variability (HRV) is presented based on the log-periodic power law 
scaling of the heart beat interval increment. Our analysis suggests mul- 
tiple DSI groups and a dynamic cascading process. A cascade model is 
presented to simulate such a property. 



1 Introduction 

The phenomenon of heart rate variability (HRV) in humans desrcibes the beat- 
to-beat, apparently random, fluctuation of the heart ratc^. HRV measured 
by the time span between ventricular contractions, known as the beat-to-beat 
RR interval (RRi), is also known to share many characteristics found in other 
natural phenomena. For example, daytime RRi in healthy humans exhibits 1/f- 
like power spectrum^, multifractal scaling'^''*, and similar increment distribution 
observed in fluid turbulence*. These characteristics may vary significantly in 
heart disease patients depending on the severity of the disease^'^. 

The origin and the generation of HRV remain the biggest challenges in the 
contemporary HRV research. Although the respiratory and vascular systems 
constantly modulate the heart rate, they do not explain the large percentage of 
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the broad-band (multifractal) signal power in HRV. For example, it is unlikely 
that this broad-band feature results directly from the output of the narrow-band 
respiratory dynamics^. Also, it is known that the level and the variability of 
blood pressure and heart rate can change significantly from upright to supine 
positions. In a 42-day long bed rest test, Fortrat et al. showed that the variation 
in blood pressure and heart before and after the test are qualitatively different, 
suggesting separate control mechanisms for generating their variability^. It is 
thus believed that a more sophisticated structure may exist, which integrates 
the feedback from receptors to create the pattern of HRV^. 

Apart from its origin, some progess on the HRV generating mechanism may 
be possible by using the discrete (lattice) multiplicative cascade model^. This 
is purely a phenomenology approach that does not prescribe to any physiology 
term. Nontheless, encouraging results were obtained that are consistent with 
the physiological data in health and in certain heart disease*. The main pur- 
pose of this work is to investigate the basis of this modeling strategy. Our 
approach is based on the scale invariant symmetry implied from the HRV 
phenomenology*'^ Since RRi cannot be defined between heart beats, it is 
appropriate to consider discrete scale invariance (DSI) in HRV. It is known that 
discrete cascade implies DSI^^""^^. Better characterization of DSI in HRV is 
thus important since it is the necessary condition for the multifractal scaling 
observed in HRV. The existence of cascade is also significant because it rep- 
resents a very different approach of the cardiovascular dynamical system from 
feedback control that is additive in principle. The idea will support the previous 
studies that a direct influence from baroreflex to multifractal HRV is unlikely^, 
as well as the need to search for a role by the higher control centers in HRV^. 

The consequence of DSI is an oscillating scaling law with a well-defined 
power law period. Such a scaling law is said to exhibit log-periodicity (LP). 
In this work, we analyzed DSI in daytime healthy HRV by searching LP in 
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the scaling of HRV. Typically, LP is "averaged out" in the process of finding 
the scaling law. Using the technique called "rephasing," this problem can be 
effectively resolved and evidence of multiple DSI groups in the healthy daytime 
RRi data was found. In light of this new result, a cascade model is constructed 
using random branching law to reproduce not only some ofthe known HRV 
plicnomcnology, but also the multiple DSI characteristics. 

The results of this work are organized in five sections. In Section 2, a brief 
review of the notion of DSI is given. The numerical procedures for identifying the 
DSI property from time series are described in Section 3. Numerical examples 
and results on daytime heart rate data sets are given in Section 4. Concluding 
remarks are given in the last Section. 

2 Discrete Scale Invariance and Rephasing 

2.1 Ideas of Discrete Scale Inveiriance in Physical Systems 

A random processes x(t) is said to possess continuous scale invariant symmetry^'^ 
if its distribution is preserved after the change of variables, t ^ Xt, x x/ fx 
where A and /Lt(A) are real numbers; i.e., 

x{t) = -x{Xt). (1) 
M 

DSI is defined when (1) only holds for a countable set of scale factors Ai, A2, ■ • •• 
Scale invariance implies power law. The power law in DSI has a log-periodic 
correction of frequency l/log(A): i.e., x{t) = F {log{t) / log{X)) where 7 = 
log/u/logA and F{x) = F{x + 1). Generally^^, one can consider x{t) = Ctf^ , 
Cf being t-dependent, and 7' = 7 + 'Inni/ log(A) is a complex number for n = 
1, 2, • • •. Novikov suggested LP in the small scale energy cascade of the turbulent 
flow^^. Sornette and co-workers showed that LP exists more generally in physical 
and financial systems, such as turbulence^^, earthquake^^, rupture^^ and stock 
market crashes^^. 
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The existence of the discrete scale factor impUes a hierarchical structure. 
This link can be simply illustrated by the middle third Cantor set with the 
scale factor A = 3. With proper rescaling, a precise copy of the set is only 
obtained with a 3-fold manification of the scale^^. If x{t) denotes the Lebesgue 
measure at scale t, the Cantor set can be modeled by (1) using A = 3 and 
jjL = X~ i°g(2)/log(3). Thus, the power law exponent of x{t) (the box dimension of 
the Cantor set^^) assumes a log-periodic oscillation of frequency 1/ log(3) about 
its mean value log(2)/ log(3). 

The hierarchical structure can be a dynamic object as a result of some time- 
dependent branching law. Such a dynamic hierarchy is believed to exist, for 
example, in the cascade paradigm of the energy exchange in fluid turbulence 
where the break-down or "branching" of large-scale vortices into ones of smaller 
scales can occur randomly in space-time with the energy re-distribution fol- 
lowing a multiplication scheme. In data analysis, the dynamic hierarchy poses 
a technical difficulty for finding the scale factor A since LP may be averaged 
out in the process of obtaining the power law. Zhou and Sornctte proposed to 
conduct averaging after rephasing or rc-aligning data segments using a central 
maximum criterion^^. Using this technique, these authors successfully extracted 
LP in turbulence and proposed the DSI symmetry and cascade. The rephasing 
technique is adopted in this work. Instead of the central maximum criterion, 
the cross-correlation property of the data segments will be used (see step (d) 
below). 

2.2 Rephasing RRi Data 

Let r(t) denote the RRi between the t^^ and {t + Vf^ heart beats. Based on 
the turbulence analogy of HRV^, we focus on the LP in the scaling exponent of 
the empirical law S{t,p) = (|Ar(T)|P) ~ r'^('^'P) where Ar(T) = Y{t + r) - r(t) 
and p is a real number. The implementation of the rephasing follows a 8-step 
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algorithm; see Fig. 1. 



(a) Divide r(t) into M nonovcrlapping segments {T^i, • • • , TZm}- 

(b) For r{t) € TZi, calculate <Si(r,p) = (|Ar(r)|P). 

(c) Apply a low-pass {K, L) Savitzky-Golay (SG) filter* to Si{T,p) and calculate 
its first derivative to obtain a r-dependent Qi{T,p) for i = 1, • • • , M. The {K, L) 
SG filter performs a K*'^ order polynomial fit over L samples^^. It can produce a 
smoothing effect in the high frequency while preserving the statistical moments 
of the signal up to the order of the filter. 

(d) Randomly select the ith segment as the base segment and compute the 
cross-correlation C^-^{n) between C,i{T,p) and Q{t + K,p) for j ^ i. 

(e) Shift the time origin of (j{T,p) by Aj, where max(C/^-^) — C-'^j^ {—Aj), so 
that the cross-correlation between Q (r) , p and the shifted Q {Aj + r, p) has a 
maximum at zero time lag. Note that Aj = for the base segment. 

(f) Average the shifted Q{t + Aj,p),j = !,•■•, M, to obtain Zk,l{t,p). 

(g) Compute the spectrum of Zk,l{t,p)- 

(h) Return to (c) with different K, L values. 

A Lomb periodogram^° ^{f) is used to estimate the spectrum of Zk,l{t,p) for 
its superiority in handling situations where noise plays a fundamental role in 
the signal, as well as its capability in handling small data set. 

Although the above algorithm provides the systematic steps to estimate the 
log-periodic component, noise in the empirical data can also generate spurious 
peaks in the Lomb periodogram. For independent Gaussian noise process, this 
problem can be analyzed by the false alarm probability^"'^^: 

Paif) = !-(!- exp(-£(/)))'" (2) 

where m is proportional to the number of points in the spectrum. The smaller 
the value Paif) is, the more likely a genuine log-periodic component exists in 
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the signal. Thus, a Lomb peak with large Pcif) suggests a chance event. Zhou 
and Sornette conducted extensive simulations and showed that (2) is in fact an 
upper bound for a number of correlated noise except for those showing long- 
term persistent correlation^^. The fractional Brownian motion (fBm) of a Hurst 
exponent greater than 0.5 is an example where (2) docs not apply. The multiple 
scaling exponents in healthy daytime HRV have been found to lie below such a 
threshold^'"^'*'*' and we will continue to use (2) in this work. 

As shown above. DSI is characterized by the frequency l/log(A) of the 
LP. However, significant Lomb peaks may only capture the higher harmonics 
k/ log(A), k ^ 1. It is therefore necessary to define the relation of the significant 
peaks. We propose a simple procedure to achieve this. First, we collect the 
significant peaks satisfying Pq < Pq for <C 1 and for different SG filter 
parameters. Second, we form a significant Lomb peak histogram (SLPH) and 
locate its local maxima. These maxima identifies the most probable frequencies 
of the log-periodic oscillation of the power law. Let such maxima be /i, • • ■ , /„. 
The last step of the procedure is to search the smallest A to minimize 



for integers fc^'s. We seek the smallest A since, with finite precision in numerical 
computing, d\ can be made arbitrarily small as A 3> 1 This minimization step 
is simple, easy to implement and, as we show below using synthetic data, it is 
also effective. 

3 Numerical Results 

3.1 DSI in Discrete Bounded Cascade 

The rephasing algorithm introduced above was first tested on synthetic data 
generated by the discrete cascade^ 
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xj{t) = nj 



[j0Jj{t) 



(4) 
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where the cascade components uij{t) are discrete-time processes given by 

ojj{t) = l+ajW (5) 

for t^^^ < t < t'j^li, fc = 1, 2, ■ • •, J = 1, ■ • • , J, and w is a zero-mean Gaussian 
random variable of variance 1. Let t^f^^ — t^^^ = Sj. The scale factor A in the 

DSI hierarchy is related to t^j''' 's by 

Sj/Sj+i = A. (6) 

To model the bounded RRi, we further assume > aj+i to assure bounded- 
ness. This model has been used in the past to simulate HRV phenomenology, 
including transition of RRi increment probability density function and multi- 
fractal scaling*. 

According to (4), we generated 30 sets of dyadic (A = 2) and triadic (A = 3) 
xj{t) with the corresponding log(aj) = (—1.6 — 0.126j) log(2) and log((7j) = 
(—1.9 — 0.093j) log(3), respectively. Each xj{t) has 8192 points and is divided 
into segments of 1024 points. Twenty-four sets of (K, L) SG filter arc defined 
based on = 3, 4. ■ • • , 7, L = 7, 9, • • • , 15. For each combination of K, L, steps 
(c) to (h) in the rephasing algorithm is repeated six times based on six different 
base segments selected in step (d) of the algorithm. This is implemented to 
avoid bias from a particular segment. Significant Lomb peaks are collected 
based on the false alarm probability Pq < 1% or C{f) > 10 and m = 256 points 
of the Lomb periodogram. The results for p = 2 is reported as no quantitative 
difference exists for p < 3. Numerical results for p > 3 show more variability 
due to poor statistics. 

FIG. 2a shows the Ci{T,p) of a particular segment of one of the dyadic x,/(t)'s. 
The log-periodic oscillation with a log-period log(2) is clearly seen. The Lomb 
periodogram of Zk,l{t,p) (step (f) above) is shown in FIG. 2b based on a 
particular choice of K, L and the dominant LP is seen to pick up the second 
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harmonics of 1/ log(2). The SLPH estimated for different SG filters over 30 sets 
of xj{t) is obtained in FIG. 3. The clustering of the local maxima at integer 
multiples of l/log(2) is evident. The minimization (3) identifies the correct 
scale factor A = 2 for the dyadic cascade. Similar results of the tradic cascade 
are also found (FIG. 3). These examples demonstrate the effectiveness of the 
proposed numerical procedures. 

3.2 DSI in Daytime Healthy HRV 

For HRV, two databases are considered. The first set (DBl) consists of 10 
ambulatory RRi recordings from healthy young adults*. These test subjects 
were allowed to conduct normal daily activities. The second set (DB2), available 
from the public domain^^, consists of 18 ambulatory RRi recordings showing 
normal sinus rhythm. The parameters used in the numerical analysis are the 
same as above except the data segment length has increased to 2048 points. The 
choice of the segment length is a balance of two factors: small segment length 
results in more segments but poorer statistics in the estimation of Ci{T,p); large 
segment length results in less segments but better estimate of Ci{T,p). We tried 
1024 points per segment and found similar results; i.e., the group averaged A 
value is similar to the ones reported in FIG. 5 below. 

The SLPH in all cases shows well positioned local maxima that can be easily 
related to the harmonics of some fundamental frequency (FIG. 4). The A values 
for DBl and DB2 are summarized in FIG. 5. It is observed that (a) there are 
non-integer scale factor A and (b) the A's cluster in the range of [3.5, 5.5] and 
the group averaged A are ~4.8 and ^4.4 for DBl and DB2, respectively. The 
noninteger A unambiguously excludes the possibility of discrete cascades with 
one scale factor. It implies more complicated branching law and multiple DSI 
groups in healthy HRV. 

Although HRV and turbulence exhibit similar phenomenology*, it is inter- 
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esting to point out the rather large A value (> 4) compared with the A ^ 2 in 
fluid turbulence^". Prom the discrete cascade viewpoint, a larger A is compati- 
ble with the "patchiness" appearance commonly observed in the RRi of healthy 
humans^'^'^'^'^° since the Wj(f)'s of the cascade will fluctuate on a longer time 
scale to create the effect. 

To model the multiple DSI in cascade HRV, the scale factor A used in (5) 
is set to be a random number so that the log-periodic oscillation of ({p) can 
vary over a range of frequencies. We generated 30 sets of xj{t) according to (4) 
using uniformly distributed A in the interval [2,6]. The simulated xj{t) exhibits 
the "patchiness" pattern observed in the RRi data (FIG. 6), and similar scaling 
characteristics found in the past^ (FIGs. 7a ~ 7c). The scaling exponent C{t,p) 
of the power law {\Axj{t)\p) ~ r^^^'^^ exhibits log-periodic oscillation that 
is captured by the well positioned local maxima in SLPH (FIGs. 7d, 7e). In 
addition, the average of the A's lies close to the group-averaged A values of DBl 
and DB2 (FIG. 5). 

4 Conclusion 

It is known that discrete cascade leads to DSI and characterized by log-periodic 
modulation of the scaling property^^'^^. Hence, the LP reported in this work 
supports the view of a cascade for the multifractal generation in HRV. It implies 
a more sophisticated process than reflex-based control mechanisms that function 
on the additive basis. It also suggests the need to search for a role by the higher 
control centers in HRV^. It is conjectured that the cascade describes the process 
which integrates the regulatory feedbacks in the cardiovascular system to create 
the pattern of HRV. 

The non-integer scale factor implies multiple DSI. This property was also 
reported in the screening competition of the growth of diffusion limited aggre- 
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gation modeP^'^^. To the best of our knowledge, this is the first instance of 
multiple DSI being reported in HRV. We do not have the better knowledge of 
its origin, except to believe it reflects the multiple time-scale control mechanisms 
in the cardiovascular dynamical system. 

It is tempting to search for the physiological correlate of the cascade, for ex- 
ample, the role of the cascade components ujj{t). Based on the spectral analysis, 
we suggested that the large time scale components ^ 1) capture mainly 

the sympatho- vagal interaction and the small time scale components (wj, j ^ 1) 
capture the parasympathetic activity"*. However, we should caution that cas- 
cade is a modeling tool derived from statistical physics. The ujj{t) can therefore 
represent the range of micro- to macroscopic processes in the cardiovascular 
dynamical system. 

A rather narrow range of the scale factor A G [3.5, 5.5] estimated from the 
two different databases implies a "stable" hierarchical structure of the cascade 
that does not vary sensitively with the details of the healthy population. The 

analysis of the identified DSI characteristics in other physiological conditions is 
currently underway and its result will be reported in the near future. 
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Figure Captions 



FIG. 1 Sketch of the numerical procedure for rcphasing. The second seg- 
ment is illustrated as the base segment and rephasing was shown for (i{t + 
■ ■ ■ ) Cm(t + Am,p) (Aj is determined at the maximum of the cross- 
correlation function between the jth and the base segments). Log-periodicity 
in Zk,l{t,p) is estimated from the Lomb periodogram. 

FIG. 2 (a) Ci{TjP) ~ {Ci{TjP)) versus log(r) taken from the synthetic dyadic 
bounded cascade. The solid line is a pure sine wave with a period of log(2) ~ 
0.693. (b) Typical Lomb periodogram of Zk,l{t, 2) (averaged over all Ci(T,p)'s). 

FIG. 3 SLPH estimated from 30 sets of (a) synthetic dyadic bounded cascade 
xj{t) and (b) triadic xj{t). The grid lines in (a) and (b) are drawn according 
to fc/log(2) and fc/log(3), A; = 1,2, • • •, respectively. 

FIG. 4 (a) SLPH of a typical data set from DBl. The local maxima /max 
are marked by (b) /max versus fc/log(4.5), k = 1,2, showing as the 

harmonics generated by the fundamental frequency l/log(4.5). The straight 
line has the slope l/log(4.5). (c) Similar to (a) based on a data set taken 
from DB2. (d) Similar to (b) based on the local maxima of (c). The straight 
line has the slope l/log(3.1). Note the local maximum between 7/log(3.1) and 
8/log(3.1) was not fitted by the harmonics of l/log(3.1). 

FIG. 5 Scale factor A's for 10 subjects in DBl, 18 subjects in DB2 and 30 sets 
of synthetic data xj{t) generated by the cascade model. The group averaged A 
values and standard deviations are superimposed and drawn as "•" and vertical 
bar, respectively. 

FIG. 6 A typical sample of xj{t) (top) and the RRi data (bottom) taken from 
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DB2. Both time series show the characteristic of "patchiness" in their fluctua- 
tion pattern. 

FIG. 7 (a) to (c) show the l//-hke power spectrum, power law, S{t,p) and the 
nonUnear ({p) of 5(t,p), respectively, of the xj{t) shown in FIG. 5; see Ref. 4 
for the similar characteristics reported for RRi data in healthy humans, (d) and 
(e) show the SLPH of two typical xj{t). Well-positioned local maxima in (d) 
and (e) capture the harmonics generated by A: -^4.4 and -^3.85, respectively. 
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This figure "figl_algorithm.jpg" is available in "jpg" format from: 
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